Spontaneous symmetry breaking of solitons trapped in a double-channel potential 
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We consider a two-dimensional (2D) nonlinear Schrodinger equation with self-focusing nonlinearity 
and a quasi-lD double-channel potential, i.e., a straightforward 2D extension of the well-known 
double-well potential. The model may be realized in terms of nonlinear optics and Bose-Einstein 
condensates. The variational approximation (VA) predicts a bifurcation breaking the symmetry 
of 2D solitons trapped in the double channel, the bifurcation being of the subcritical type. The 
predictions of the VA are confirmed by numerical simulations. The work presents the first example 
of the spontaneous symmetry breaking (SSB) of 2D solitons in any dual-core system. 
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I. INTRODUCTION 

The creation of Bose-Einstein condensates (BECs) in 
vapors of alkali metals has opened up a unique opportu- 
nity to investigate nonlinear interactions of atomic mat- 
ter waves. In particular, an important physical prob- 
lem is to develop methods allowing one to create and 
control matter- wave solitons in the experiment. One- 
dimensional (ID) dark [l[, bright Q, and gap- mode Q 
solitons have already been created. 

Another fascinating phenomenon, viz., the Josephson 
effect in BEC, was observed in condensates loaded in a 
double- well potential. Tunneling through a potential bar- 
rier usually occurs in quantum systems on a nanoscopic 
scale, while the Josephson effect features the tunneling 
of macroscopic wave functions describing intrinsically co- 
herent states [i| . This phenomenon has been observed in 
sundry systems, such as a pair of superconductors sep- 
arated by a thin insulator (the Josephson effect proper) 
Q , and two reservoirs of superfluid helium connected by 
nanoscopic apertures [1, 0] • Recently, the first successful 
implementation of a bosonic Josephson junction formed 
by two weakly coupled BECs in a macroscopic double- 
well potential was reported [8(. In contrast to hitherto 
realized Josephson systems in superconductors and su- 
perfluids, interactions between tunneling particles plays 
a crucial role in the junction implemented in the BEC 
setting, the effective nonlinearity induced by the inter- 
actions giving rise to new dynamical regimes in the tun- 
neling. In particular, anharmonic Josephson oscillations 
were predicted [j| |T(| [H| , provided that the initial pop- 
ulation imbalance in the two potential wells falls below 
a critical value. The dynamics change drastically for the 
initial population difference exceeding the threshold of 
the macroscopic quantum self-trapping and thus inhibit- 
ing large-amplitude Josephson oscillations Q, [H, [HI ]. 
These two different regimes have been experimentally 
demonstrated in the BEC-based Josephson-junction ar- 
rays Q [HQ- 

This dynamics can be well explained by means of a 



simple model derived from the Gross-Pitaevskii equation 
(GPE). Two equations for the self- interacting BEC am- 
plitudes, linearly coupled by tunneling terms, describe 
the dynamics in terms of the inter-well phase difference 
and population imbalance. As mentioned above, the 
nonlinearity specific to BEC gives rise to the "macro- 
scopic quantum self-trapping" effect, in the form of a 
self- maintained population imbalance Q, [TH, [ill [2(| ■ In 
order to derive a reduced two-state model, one needs to 
find eigenstates of the corresponding GPE and perform 
stability analysis for them. Such analysis can be readily 
performed for the square-shaped double-w ell p otential, 
where analytic solutions are available [HI, l2Ct [2l| . In 
this case, the point of the symmetry-breaking bifurca- 
tion, where asymmetric solutions emerge, can be found 
exactly. 

The present paper addresses the symmetry-breaking 
bifurcation and the existence and stability of asymmetric 
states in a two-dimensional (2D) system, which is a direct 
extension of the familiar double- well model QQ. The 
corresponding potential is shown in Fig. [T] below. It fea- 
tures the double-square- well shape in the x direction, be- 
ing uniform along y. To our knowledge, the spontaneous 
symmetry breaking (SSB) has never been studied in 2D 
systems before. By means of the variational approxima- 
tion (VA), we will find regions where stable asymmetric 
solitons exist, and the prediction will be then confirmed 
by direct simulations. 

The paper is organized as follows. The model in in- 
troduced in Sec. II, where its physical interpretations 
are outlined too, in terms of BEC and nonlinear optics. 
Then, in Sec. Ill, we derive variational equations and 
analyze their solutions, which predict the SSB of a sub- 
critical type. At the end of that Section, we compare the 
result with those for the CW (continuous- wave, i.e., as 
a matter of fact, one-dimensional) states, for which the 
SSB bifurcation is of a different type, being supercriti- 
cal. In Sec. IV, we compare predictions of the VA with 
numerical results, and Sec. V concludes the paper. 
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FIG. 1: The shape of the quasi-one-dimensional double-well 
potential, U(x,y). 



II. THE MODEL 

The starting point is the 2D equation in a normalized 
form, with the self-attracting cubic nonlinearity, 

m + \ (* M + *„,) - u(x)* + lif* = o, (i) 

where the quasi-lD double-well potential is taken as 



U(x) 



0, |at| < \L and \x\ > \L 
\L < \x\ <\L + D," 



D. 



-U , 



(2) 



with D, Uq and L being, respectively, the width and 
depth of each well, and the width of the barrier between 
them, see Fig. [TJ 

Equation ([1} admits three different physical interpre- 
tations. In terms of BEC, it is the Gross-Pitaevskii 
equation (GPE), for a pancake-shaped (nearly flat) self- 
attractive condensate (in fact, this implies 7 Li , in which 
the scattering length can be readily made slightly nega- 
tive by means of the Feshbach resonance Q), with po- 
tential @ acting in the (x,y) plane. In nonlinear optics, 
the evolutional variable, t, is actually the propagation 
distance (usually denoted by z). Then, one possible in- 
terpretation is that Eq. ([TJ) is the nonlinear Schrodinger 
equation (NLSE) which, adopting the ordinary paraxial 
approximation for the diffraction, governs the stationary 
transmission of an optical signal in the spatial domain 
(i.e., in the bulk medium) with the self- focusing Kerr 
nonlinearity; then, the double-channel potential corre- 
sponds to two waveguiding slabs embedded in the 3D 
medium. An alternative optical interpretation is valid 
in the temporal domain, where Eq. ([T]) describes the 
light propagation in a nonlinear planar waveguide with 
a pair of embedded guiding channels. In the latter case, 
x is the transverse coordinate, while y is actually the re- 
duced time, t — z/V (t is physical time, and V is the 
mean group velocity of the carrier wave), assuming that 



the group-velocity dispersion in the planar waveguide is 
anomalous. 

According to the above interpretations, possible local- 
ized solutions of equation JTJ), if it is considered as the 
GPE, will be interpreted as matter-wave solitons. If the 
equation is realized as the NLSE in optical media, the lo- 
calized solutions will be either spatial or spatiotemporal 
solitons (in the bulk and planar waveguide, respectively; 
a review of the topic of multidimensional optical solitons 
can be found in Ref. [22j). 

We are looking for stationary solutions as ^(x,y,t) = 
e~ l ^<&{x, y), where real function &(x,y) satisfies equa- 
tion 



1 , 
2 



$„„) - U(x)$ + $ 3 = 0, 



(3) 



which can be derived from the Lagrangian, 



Is 



dxdy 



-U(x)<f> 2 



&i + <F:, 



(4) 



Our variational approximation (VA) will be based on this 
Lagrangian. 



III. VARIATIONAL APPROXIMATION 

A. Symmetry breaking of the two-dimensional 
solitons 

To apply the VA (a detailed account of the technique 
can be found in review [23}), we assume a situation with 
two very narrow and deep channels, separated by a broad 
barrier, i.e., D -C L in Eq. ([2]). The ansatz describing 
the soliton field configuration consists of two parts. First, 
inside each channel, i.e., in regions |ar =p (L + D) /2| < 
D/2, we adopt the trial function 



&(x, y) = A± cos I 7r 



xt (L + D)/2 
D 



exp 



y 



2W 2 



where A± and W are three real variational parameters. 
In this expression, we assume that the wave function has 
different amplitudes but equal longitudinal widths, W, 
in both channels. In direction x, ansatz ([5]) emulates the 
ground state of a quantum-mechanical particle in an in- 
finitely deep potential box, therefore it vanishes at edges 
of the channel. In direction y, the ansatz is assumed 
to be a self-trapped soliton, approximated by the Gaus- 
sian. The form of the ansatz outside the channels is also 
suggested by quantum mechanics, emulating a superpo- 
sition of exponentially decaying ground-state wave func- 
tions behind the edges of deep potential boxes, 



$(x,y) = J2 A ± ex P ( -V-%i* 
+,- ^ 



XT 



L + D 



y 



2W 2 



(6) 
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where amplitudes A± and width W are the same as in 
Eq. ©. 

The substitution of the inner and outer parts of the 
ansatz, Eqs. ([5|) and (J6j> , into Lagrangian (j4]), and sub- 
sequent integration over x and y lead to the following 
effective Lagrangian: 



hence parameters N± are proportional to the populations 
in the channels, while v determines the population im- 
balance. In terms of this notation, effective Lagrangian 
([7]) transforms into 



-L 



off 



E 



8W 



(7) 



Here, we have adopted the Thomas-Fermi approxima- 
tion in the x direction (but not along y), by omitting 
the corresponding kinetic-energy term, — (l/2)$> 2 , in the 
Lagrangian density. 

Note that ansatz ([5]) includes the cosine trial function 
with a constant width. This assumption is relevant for 
a deeply bound quantum state (as mentioned above, the 
ansatz was modeled on the pattern of the ground state 
in the infinite deep box) , but not when the energy eigen- 
value \fj,\ is small. Indeed, if one tries to predict, by 
means of this ansatz, a bound state of a particle in a 
finite-depth rectangular potential well in ordinary (lin- 
ear) quantum mechanics, one would arrive at a conclu- 
sion that the bound state appears only starting from a 
minimum finite value of Uq, which is, as a matter of fact, 
the kinetic energy in the x direction. A commonly known 
exact result is that there is no threshold for the existence 
of a bound state in any symmetric potential well, even if 
it is arbitrarily shallow. However, this unphysical effect 
disappears in the Thomas-Fermi approximation, which 
was adopted above. 

To simplify the notation, we introduce new parameters: 

e = (i + Uq (8) 
(notice that e may be both positive and negative), 



A = (2/D) v/^exp (-y/=2jl(L + D) 



and N± = (3/4 V^) A\W . Additionally, we define 



N 



N 4 



A_ 



4VA 



A+ - A_ 
4VA 



(9) 



(10) 



Norms of the wave function (which are proportional to 
the numbers of trapped atoms) in the two channels are, 
according to ansatz ©, 



f + oo 



dy 



±{D+L/2) 



±L/2 



dx{<S>{x,y))' 



= (V5F/2) A 2 ± DW, 



4VA f-L V 2 8W2 m J 2 



e A A y/\ A 2 + , 



2 8W 2 2 W 



+ \\jN 2 -v 2 . (11) 



This Lagrangian gives rise to variational equations 
dL/dW — dL/dN — dL/dv — 0, i.e., respectively, 



N 



2VA(A 2 + v 2 ) 
1 1 V\N XN 
2 £ ~ 8W 2 + W + y/N 2 - v 2 

( Va a \ 



W y/N 2 - V 2 



- W, (12) 
= 0, (13) 
= 0. (14) 



Equation (jl4|) has two solutions: v = 0, which, accord- 
ing to Eq. (|10p corresponds to symmetric solitons, and 
asymmetric ones, with v ^ 0. 

In Eqs. (|12[) - (fT4"|) , it is relevant to consider N (which 
is proportional to the total number of atoms) as a given 
parameter, and e [which is, as the matter of fact, the 
chemical potential, see Eq. (J8J)] as an unknown. In this 
manner, for each value of N we can find the correspond- 
ing vales of A, e, v and W [system (fl2|) - (fT4|) contains 
only three equations, but, as seen from definitions ([8]) 
and ((9J, e and A are not independent]. 

Being interested in asymmetric solutions, with v ^ 0, 
and substituting expression (TT2"|) for W in Eq. (TLT)) , we 
arrive at an equation which determines the population 
imbalance, v 1 as a function of N: 



2\/N 2 - v 2 (N 2 + v 2 ) = N. 



(15) 



Notice that this cubic equation for v 2 does not contain 
A, W or e. The most essential issue is when the spon- 
taneous symmetry breaking (SSB) occurs, i.e., at what 
value of a nontrivial solution for v appears in Eq. (| 1 5|) . 

Straightforward analysis of Eq. (|T5)) demonstrates 
that, with the increase of A, a pair of physical (real) 
solutions for v appears through a tangent (alias saddle- 
center) bifurcation at A 



A min = (1/2)^3^3/2 a 
0.958. At a slightly larger critical value, A max = 1, a 
subcritical pitchfork bifurcation takes place, giving rise 
to solutions splitting off from v = 0. The entire bi- 
furcation diagram for Eq. (|15p is displayed in Fig. 
According to general principles of the bifurcation the- 
ory [24j, the picture implies that the symmetric solution, 
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FIG. 2: (Color online). The subcritical symmetry-breaking 
bifurcation in the double-channel model, as predicted by 
the variational approximation through Eq. (|15[) . Unstable 
branches of the solutions are shown by dashed curves. 



v = 0, is stable in interval < N < 1, and unstable 
for N > 1. Simultaneously, branches of the asymmetric 
solutions (those with u ^ 0) are unstable as long as they 
go backward, and become stable after they turn forward 
at point N = N ni i n . 

This subcritical bifurcation is qualitatively similar to 
that explored earlier the model of dual-core nonlinear op- 
tical fibers, which is based on a pair of linearly coupled 
one-dimensional NLSEs for amplitudes of the electromag- 
netic waves in the two cores [23|, [25| . On the other hand, 
in an apparently more complex model of parallel-coupled 
fiber Bragg gratings, that amounts to a system of four 
equations with linear and nonlinear couplings, the SSB 
is simpler, being supercritical (the emerging branches of 
asymmetric solutions immediately go forward and are 
stable everywhere) (26|. 



B. Comparison with the one-dimensional 
(continuous-wave) case 

It is relevant to compare the above results for the SSB 
of solitons in the 2D model to what can be predicted in 
the ID counterpart of the model by the same type of the 
VA. The latter corresponds to the ansatz based on Eqs. 
(0 and ©, but with W — oo (in other words, this is a 
CW state in terms of the 2D model). Then, Lagrangian 
J7]) reduces to 

L cS = const ■ ]T {^eA 2 ± + + 2XA + A^j . 

(16) 

Straightforward manipulations with the variational equa- 
tions generated by this Lagrangian, dL c g/dA+ = 
dL e g/dA_ — 0, yield a final relation for asymmetric 
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N 

FIG. 3: (Color online). The supercritical symmetry-breaking 
bifurcation for the continuous-wave (CW) states, in the ID 
model. 



states: 

(A 2 + - A 2 _f = (A\ + A 2 _) 2 - 2 (16A/3) 2 . (17) 

In this context, N = (3/32A) (A\ + A 2 _) is again pro- 
portional to the norm, and v = (3/32X)(A 2 i _ — A 2 _) may 
be considered as a measure of the asymmetry. The pur- 
port of Eq. (fl7|) is that it predicts a critical value of N 
at which the asymmetric solutions emerge, N CI = 1. As 
shown in Fig. [3J a principal difference of the SSB for the 
CW states from its counterpart for the solitons (see Fig. 
[2]) is that the present bifurcation is a supercritical one. 

It is relevant to mention that the SSB for the CW so- 
lutions of the dual-core-fiber model is also supercritical 
[13], on the contrary to the subcritical bifurcation for 
solitons in the same model [25| (the bifurcation of the 
CW solutions may become subcritical if the nonlinear- 
ity is saturable, rather than cubic Finally, it is 

also necessary to mention that all the CW states, con- 
sidered as quasi- ID solutions of the 2D model with the 
self- attracting nonlinearity, are unstable against modu- 
lational perturbations (while solitons may be completely 
stable in long simulations, see below). 



IV. NUMERICAL RESULTS 

To verify the above predictions, we solved Eq. §3§ nu- 
merically, using the imaginary-time relaxation method 
with a fourth-order Runge-Kutta algorithm. The accu- 
racy of the numerical code was tested by varying com- 
putational parameters, namely the mesh density, domain 
size, and time step. These parameters were then fixed at 
values for which further increase of the accuracy would 
not lead to a visible change in the final results [28| . This 
procedure was applied every time when the physical pa- 
rameters L, D, N were varied. 



5 



The stability of solitons produced by this method was 
then tested by direct integrations of perturbed states in 
real time. The perturbation was introduced, multiply- 
ing the wave function by a symmetry-breaking factor, 
^ — > ^(1 + ay/L). Perturbations with a = 0.05, which 
are actually large, were not able to destroy solitons that 
were identified as stable ones. On the contrary, much 
small perturbations (for instance, with a = 0.002) were 
sufficient to demonstrate instability of those solutions 
which are unstable, after propagation time t = 100. 

In this work, we did not attempt to identify stability 
regions for the solitons by computing their eigenvalues 
in terms of linearized equations for small perturbations, 
so, in this sense, the stability borders are not completely 
rigorous ones. Nevertheless, the distinction between un- 
stable and stable solitons revealed by the simulations is 
very clear, and, on the other hand, the identification of 
the stability by dint of direct simulations corresponds 
to experimental conditions, where solitons are subject to 
various perturbations of a finite size. 

The numerical results are summarized in Figs. 0J[7J 
Figure Q] displays a typical dependence of the global 
asymmetry coefficient for the numerically found soliton 
solutions, defined as 



77,-1- — 71- 

n 



/-to <& 2 dxdy 



(18) 



on the total norm, n. The way the figure was gener- 
ated (through direct simulations converging to stationary 
states) made is possible to display only stable branches of 
the solutions, both symmetric and asymmetric ones. Al- 
though the unstable branches are missing, there is little 
doubt that the full SSB diagram corresponding to Fig. 0] 
is of the generic subcritical type. In particular, a bistabil- 
ity (hysteretic) region, where symmetric and asymmetric 
solitons coexist and are simultaneously stable, is obvious 
in the figure. Thus, the picture suggested by the numer- 
ical results is fully consistent with the predictions of the 
VA shown in Fig. [21 

An example of coexisting symmetric and asymmetric 
solitons is shown in Fig. for the value of the norm 
n = 4.05. As both solitons contain equal total numbers 
of atoms, the asymmetric one has smaller width in both 
x and y directions, as its shape provides for stronger ef- 
fective self-attraction. On the the other hand, the shape 
of the symmetric soliton features a fair amount of tun- 
neling between the channels. Figures a) and b) were 
generated by direct numerical simulations, and Figs. [5] c) 
and d) were obtained from the VA. It is observed that 
the agreement between the numerical and variational re- 
sults is quite good. For a more detailed comparison, we 
display cross-sections of both solitons in panels e) and f). 

Figure [5] illustrates the stability of various solitons in 
numerical simulations. Below the bifurcation point (for 
n = 3.9), we present stable evolution of the symmetric 
state, and above the bifurcation we demonstrate a stable 
asymmetric state for n = 4.15. Also, for norm n = 4.15, 




FIG. 4: (Color online). Imbalance in the norm between half- 
planes x > and x < 0, defined as per Eq. (|18[) for numeri- 
cally found stationary soliton solutions, versus the total norm. 
The continuous and dashed lined lines show, respectively, nu- 
merically found stable symmetric and asymmetric solutions 
(unstable solutions were not generated by the numerical pro- 
cedure). Parameters are: L = 1, D = 1, and Uo = 1. 



we display the evolution of the unstable symmetric state. 
It is worthy to note that this unstable state does not sim- 
ply relax to the stable one, but rather performs persistent 
oscillations between symmetric and asymmetric shapes. 

In Fig. [7J we present comparison of the analytical re- 
sults, obtained by means of the VA for the family of 
asymmetric solitons, with numerical findings. The solid 
line is the v{N) dependence corresponding to the stable 
upper branch of the plot in Fig. [21 whereas the other 
two lines show two sets of the corresponding numerical 
results, generated as said in the caption to Fig. [7J Gener- 
ally, the solid line (variational solution) is shifted towards 
smaller values of N. 



V. CONCLUSIONS 

In this work, we have introduced a physical model that 
gives rise to the first example of the SSB (spontaneous 
symmetry breaking) of 2D solitons in a dual-core system; 
previously, this effect was studied in detail, but solely in 
ID settings. Our model is based on the 2D nonlinear 
Schrodinger equation with the self-focusing cubic non- 
linearity and a quasi-lD double-channel potential. The 
model applies to the description of spatial optical soli- 
tons in a bulk medium with two waveguiding slabs em- 
bedded into it, or spatiotemporal solitons in a planar 
waveguide into which two guiding channels were inserted. 
The same model may also be interpreted as the Gross- 
Pitaevskii equation for a pancake-shaped Bose-Einstein 
condensate trapped around two attractive parallel light 
sheets. By means of the VA (variational approximation), 
we have predicted the SSB bifurcation for the 2D soli- 
tons supported by the double channel. The bifurcation 
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FIG. 5: (Color online). Examples of numerically found asym- 
metric (a) and symmetric (b) soliton solutions of Eq. ((3]), for 
n = 4.05. Other parameters are the same as in Fig. [4] The 
gray shading indicates the position of the two potential wells. 
Panels (c) and (d) show the corresponding solutions as pre- 
dicted by the variational approximation. In panels (e) and (f ) , 
the cross-sections of the solutions along y — are compared: 
The dashed and solid curves correspond to the numerical and 
variational solutions, respectively. 




FIG. 7: (Color online). Comparison of the variational pre- 
diction for the asymmetric solitons (solid line) with respec- 
tive numerical results (dotted and dashed lines). The dashed 
curve corresponds to the same parameters as in Fig. |4j while 
the dotted one is generated by the numerical solution of Eq. 
([3]) with the distance between the wells increased to L = 2. 



was predicted to be subcritical (unlike its counterpart in 
the continuous- wave ID model). The predictions of the 
VA are well corroborated by numerical solutions. 




Wit b) 




FIG. 6: Typical examples of the evolution of perturbed sta- 
tionary states produced by numerical simulations of Eq. {1} . 
(a): A stable symmetric state for n = 3.9. (b) and (c): Sta- 
ble asymmetric and unstable symmetric states for n = 4.15. 
Other parameters are as in Fig. [4] The evolution time is 
t = 300. 
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